In [1]:
# necessary packages #

#using Pkg
#Pkg.add("Distances")
using Distributions
using Random
using Distances
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using ProgressMeter
using JLD2
In [2]:
include("../../util.j")
Out[2]:
colnorm (generic function with 1 method)
In [3]:
# unnecessary packages #

#using Pkg
#Pkg.add("UnicodePlots")
using UnicodePlots   # check the structure of the sparse matrix
using BenchmarkTools

using StatsPlots
using MCMCChains
using PrettyTables
In [4]:
#using Pkg
#Pkg.add("ProgressMeter");
In [5]:
@load "../data/sim3data.jld";
In [6]:
# priors #
q = 10; p = 3; K = 10
μβ = spzeros(p, q); inv_Vr = spzeros(p, p);
μΛ = spzeros(K, q); inv_VΛ = spzeros(K, K)
aΣ = 2; bΣ = fill(1.0, q);
ϕU = 300 / sqrt(2); ϕL = 3 / sqrt(2);
inv_Lr = spzeros(p, p); inv_LΛ = spzeros(K, K);
In [7]:
# Some data preparations #

m = 10; n = length(S);                                        # number of nearest neighbor                       
NN = BuildNN(coords_ord[:, S], m, 1.0);                            # build nearest neighbor 
nnIndx_col = vcat(NN.nnIndx, 1:n);                            # the index of columns
nnIndx_row = zeros(Int64, 0);                                               
for i in 2:m
    nnIndx_row = vcat(nnIndx_row, fill(i, i-1))
end
nnIndx_row = vcat(nnIndx_row, repeat((m + 1):n, inner = m), 1:n);  # the index of rows
In [15]:
# preallocation #
μ_m = [Array{Float64, 2}(undef, length(M_ind[i]), q) for i in 1:q];
nIndx = length(NN.nnIndx);
A = [Array{Float64}(undef, nIndx) for i in 1:K];
D = [Array{Float64}(undef, n) for i in 1:K];
I_A = [spzeros(n, n) for i in 1:K];
A_new = [Array{Float64}(undef, nIndx) for i in 1:K];
D_new = [Array{Float64}(undef, n) for i in 1:K];
I_A_new = [spzeros(n, n) for i in 1:K];
Ystar = vcat(Y_ord[S, :], inv_Lr * μβ, inv_LΛ * μΛ);             # will be updated after imputing missing response
Xstar = vcat([X_ord[S, :] spzeros(n, K)], [inv_Lr spzeros(p, K)], 
    [spzeros(K, p) inv_LΛ]);      
bstar = fill(0.0, q); astar = aΣ + 0.5 * n;
μγstar = vcat(μβ, μΛ); invVγstar = fill(0.0, p + K, p + K);
Y_Xm = spzeros(n + p + K, q);

MCMC sampling algorithm Q1: priors for $\nu_i$ Q2: $\phi_i$ may not be consistant, since the order can change

In [16]:
# Preallocation for MCMC samples and Initalization #
N_sam = 20000;
γ_sam = Array{Float64, 3}(undef, p + K, q, N_sam + 1);
Σ_sam = Array{Float64, 2}(undef, q, N_sam + 1);
F_sam = Array{Float64, 3}(undef, n, K, N_sam);
Y_m_sam =  [Array{Float64, 2}(undef, length(M_ind[i]), N_sam) for i in 1:q];
A_sam = Array{Float64, 1}(undef, N_sam); # acceptance rate
lh_old = 1; lh_new = 1;     # record the likelihood for updating ranges

ϕ_sam = Array{Float64, 2}(undef, K, N_sam + 1);

γ_sam[:, :, 1] = fill(0.0, p + K, q);
Σ_sam[:, 1] = fill(0.2, q);
ϕ_sam[:, 1] = fill(6.0, K);

precond_D = Array{Float64, 1}(undef, K * n);
inv_sqrt_Σ_diag = Array{Float64, 1}(undef, q);
RWM_scale = 0.3;                                              # random-walk metropolis step size scale
In [17]:
# for loop for MCMC chain #
Random.seed!(123);
prog = Progress(N_sam, 1, "Computing initial pass...", 50)
for l in 1:N_sam
    # Build the matrix D_Sigma_o^{1/2} #
    inv_sqrt_Σ_diag = 1 ./ (sqrt.(Σ_sam[:, l]));
    invD = Diagonal(vcat([inv_sqrt_Σ_diag[index_S_M[i, :] .> 0] for i in 1:N if (index_S[i] > 0)]...));
    if l == 1
        for k in 1:K
            getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, ϕ_sam[k, l], 0.5, A[k], D[k]);
            I_A[k] = sparse(nnIndx_row, nnIndx_col, vcat(-A[k], ones(n)));
        end
    end
    Ytilde = vcat(invD * vcat([Y_ord[S_ind[ind], ind] - X_ord[S_ind[ind], :] * γ_sam[1:p, ind, l] 
                                for ind in 1:q]...)[perm_ind], zeros(K * n));
    Xtilde = vcat(invD * kron(sparse(transpose(γ_sam[(p + 1):(p + K), :, l])), 
                            sparse(1:N, 1:N, ones(N)))[obs_ind, 
                            vcat([S .+ (ind - 1) * N for ind in 1:K]...)][perm_ind, :],
             blockdiag([Diagonal(1 ./ sqrt.(D[ind])) * I_A[ind] for ind in 1:K]...));
                  
    # use LSMR to generate sample of F #       
    nsam = length(Ytilde);
    Precond_D = colnorm(Xtilde);
    F_sam[:, :, l] = reshape(Diagonal(1 ./ Precond_D) * lsmr(Xtilde * Diagonal(1 ./ Precond_D), 
                            Ytilde + rand(Normal(), nsam)), :, K);               
                    
    # impute missing response  over S#
    Xstar[1:n, (p + 1):(p + K)] = F_sam[:, :, l];        # update matrix Xstar with F
    for ind in 1:q
        #mul!(μ_m[ind], Xstar[M_Sind[ind], :], γ_sam[:, :, l]); #Xstar[M_Sind[ind], :] * γ_sam[:, ind, l]
        Y_m_sam[ind][:, l] = Xstar[M_Sind[ind], :] * γ_sam[:, ind, l]+    # μ_m[ind][:, ind] + 
            rand(Normal(0.0, sqrt(Σ_sam[ind, l])), length(M_ind[ind]));
    end
                        
    # use MNIW to sample γ Σ #
    for ind in 1:q
        Ystar[M_Sind[ind], ind] = Y_m_sam[ind][:, l]   # update Ystar with imputed response
    end
    
    invVγstar = cholesky(Xstar'Xstar);
    mul!(μγstar, transpose(Xstar), Ystar); μγstar = invVγstar.U \ (invVγstar.L \ μγstar);
    Y_Xm = Ystar - Xstar * μγstar;      # maybe improve?
    bstar = [bΣ[ind] + 0.5 * (norm(Y_Xm[:, ind])^2) for ind in 1:q];
    Σ_sam[:, l + 1] = [rand(InverseGamma(astar, bstar[ind]), 1)[1] for ind in 1:q];          # sample Σ
    γ_sam[:, :, l + 1] = (invVγstar.U \ reshape(rand(Normal(), (p + K) * q), (p + K), q)) * 
                    Diagonal(sqrt.(Σ_sam[:, l + 1])) + μγstar;          # sample γ    
    # use metropolis-hasting to update range
    ϕ_sam[:, l + 1] = ϕ_sam[:, l] + RWM_scale * rand(Normal(), K); # propose next sample point
                    
    
    if all(x -> (x < ϕU && x > ϕL), ϕ_sam[:, l + 1])
        lh_old = - 0.5 * sum([sum(log.(D[k])) + 
                              norm((I_A[k] * F_sam[:, k, l]) ./ sqrt.(D[k]))^2 for k in 1:K]);
        for k in 1:K
            getAD(coords_ord[:, S], NN.nnIndx, NN.nnDist, NN.nnIndxLU, 
                  ϕ_sam[k, l + 1], 0.5, A_new[k], D_new[k]);
            I_A_new[k] = sparse(nnIndx_row, nnIndx_col, vcat(-A_new[k], ones(n)));
        end
        lh_new = - 0.5 * sum([sum(log.(D_new[k])) + 
                              norm((I_A_new[k] * F_sam[:, k, l]) ./ sqrt.(D_new[k]))^2 for k in 1:K]);
        A_sam[l] = min(exp(lh_new - lh_old), 1.0);
        if rand(1)[1] < A_sam[l]
            I_A = copy(I_A_new); D = copy(D_new);
        else
            ϕ_sam[:, l + 1] = ϕ_sam[:, l]; # Don't update
        end         
    else
        A_sam[l] = 0.0;
        ϕ_sam[:, l + 1] = ϕ_sam[:, l];   # Don't update when falling out of the supports
    end                      
    
    next!(prog) # monitor the progress
end

Computing initial pass...100%|██████████████████████████████████████████████████| Time: 0:54:018:40
In [18]:
mean(A_sam)
Out[18]:
0.21475224500961534

Posterior prediction

In [19]:
# prediction preparison
#N_pre_burn = Integer(trunc(0.75 * N_sam));
#M_ind = setdiff(1:N, S); NM = length(M_ind)
#F_M_sam = Array{Float64, 3}(undef, NM, K, N_sam - N_pre_burn + 1);
#Y_M_sam = Array{Float64, 3}(undef, NM, q, N_sam - N_pre_burn + 1);

# construct Atilde Dtilde #

#using RCall
#@rput coords_ord
#@rput S
#@rput m
#R"""
#library("RANN")
#nn_mod_ho <- nn2(t(coords_ord[, S]), t(coords_ord[, -S]), k = m)
#"""
#@rget nn_mod_ho
#Atilde = Array{Float64}(undef, NM * m); Dtilde = Array{Float64}(undef, NM);
#MnnIndxLU = collect(1:m:(NM * m + 1));
#MnnIndx = vec(nn_mod_ho[:nn_idx]');
In [20]:
#for i in N_pre_burn:N_sam
#    for j in 1:K
#        # update F
#        getAD(coords_ord[:, S], MnnIndx, vec(nn_mod_ho[:nn_dists]'), MnnIndxLU, 
#             ϕ_sam[j, i + 1], 0.5, Atilde, Dtilde)
#        AtildeM = sparse(repeat(1:NM, inner = m), MnnIndx, Atilde, NM, n);
#        F_M_sam[:, j, (i - N_pre_burn + 1)] = AtildeM * F_sam[:, j, i] + sqrt.(Dtilde) .* rand(Normal(), NM)
#    end 
#    # update Y
#    Y_M_sam[:, :, (i - N_pre_burn + 1)] = X_ord[M_ind, :] * γ_sam[1:p, :, i + 1] + 
#        F_M_sam[:, :, (i - N_pre_burn + 1)] *  γ_sam[(p + 1):(p + K), :, i + 1] + 
#        transpose(rand(MvNormal(Σ_sam[:, :, i + 1]), NM))
#end

MCMC Chain check

In [21]:
β_pos_sam = Array{Float64, 3}(undef, N_sam + 1, p * q, 1);
β_pos_sam[:, :, 1] = hcat([γ_sam[i, j, :] for i in 1:p, j in 1:q]...);
β_chain = Chains(β_pos_sam);
pβ = plot(β_chain)
Out[21]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -4 -3 -2 -1 0 1 2 Param1 Iteration Sample value -4 -2 0 2 0.0 0.1 0.2 0.3 0.4 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -6 -5 -4 -3 -2 -1 0 Param2 Iteration Sample value -6 -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 8 Param3 Iteration Sample value 0 2 4 6 8 0.0 0.5 1.0 1.5 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -8 -6 -4 -2 0 Param4 Iteration Sample value -8 -6 -4 -2 0 0.0 0.1 0.2 0.3 Param4 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 Param5 Iteration Sample value 0 1 2 3 0.0 0.5 1.0 1.5 Param5 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 Param6 Iteration Sample value 0 2 4 6 0.0 0.5 1.0 1.5 Param6 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -1 0 1 2 Param7 Iteration Sample value -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 0.5 Param7 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 Param8 Iteration Sample value 0 1 2 3 0.0 0.5 1.0 1.5 Param8 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -12 -10 -8 -6 -4 -2 0 Param9 Iteration Sample value -12 -10 -8 -6 -4 -2 0 0.0 0.5 1.0 1.5 Param9 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -9 -6 -3 0 3 Param10 Iteration Sample value -10 -5 0 5 0.00 0.05 0.10 0.15 Param10 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 Param11 Iteration Sample value -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Param11 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -1.0 -0.5 0.0 0.5 1.0 Param12 Iteration Sample value -1.0 -0.5 0.0 0.5 1.0 1.5 0.00 0.25 0.50 0.75 1.00 1.25 Param12 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -4 -3 -2 -1 0 1 2 Param13 Iteration Sample value -4 -2 0 2 0.00 0.05 0.10 0.15 0.20 0.25 Param13 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -6 -5 -4 -3 -2 -1 0 Param14 Iteration Sample value -6 -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param14 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -4 -3 -2 -1 0 Param15 Iteration Sample value -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param15 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -4 -2 0 2 4 Param16 Iteration Sample value -4 -2 0 2 4 0.0 0.1 0.2 0.3 Param16 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 5 Param17 Iteration Sample value 0 1 2 3 4 5 0.00 0.25 0.50 0.75 1.00 1.25 Param17 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 8 Param18 Iteration Sample value 0 2 4 6 8 0.00 0.25 0.50 0.75 1.00 1.25 Param18 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -2 0 2 4 6 Param19 Iteration Sample value -2.5 0.0 2.5 5.0 7.5 0.00 0.05 0.10 0.15 0.20 Param19 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 5 6 Param20 Iteration Sample value 0 1 2 3 4 5 6 0.00 0.25 0.50 0.75 1.00 Param20 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -10.0 -7.5 -5.0 -2.5 0.0 Param21 Iteration Sample value -10.0 -7.5 -5.0 -2.5 0.0 0.00 0.25 0.50 0.75 1.00 Param21 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -8 -6 -4 -2 0 Param22 Iteration Sample value -8 -6 -4 -2 0 0.0 0.1 0.2 0.3 Param22 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -3 -2 -1 0 Param23 Iteration Sample value -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param23 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 1 2 3 4 Param24 Iteration Sample value 0 1 2 3 4 5 0.0 0.5 1.0 1.5 Param24 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -4 -2 0 2 Param25 Iteration Sample value -6 -4 -2 0 2 0.0 0.1 0.2 0.3 Param25 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 Param26 Iteration Sample value 0 2 4 6 0.00 0.25 0.50 0.75 1.00 1.25 Param26 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 Param27 Iteration Sample value 0 2 4 6 8 0.0 0.5 1.0 1.5 Param27 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -1 0 1 2 3 4 Param28 Iteration Sample value -2 0 2 4 0.0 0.1 0.2 0.3 0.4 0.5 Param28 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -5 -4 -3 -2 -1 0 Param29 Iteration Sample value -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param29 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -6 -5 -4 -3 -2 -1 0 Param30 Iteration Sample value -6 -5 -4 -3 -2 -1 0 0.0 0.5 1.0 1.5 Param30 Sample value Density
In [22]:
β
Out[22]:
3×10 Array{Float64,2}:
  1.0  -1.0    1.0  -0.5   2.0  -1.5   0.5   0.3  -2.0   1.5
 -5.0   2.0    3.0  -2.0  -6.0   4.0   5.0  -3.0   6.0  -4.0
  8.0   6.9  -12.0   0.0  -4.0   7.7  -8.8   3.3   6.6  -5.5
In [23]:
Λ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, K * q, 1);
Λ_pos_sam[:, :, 1] = hcat([γ_sam[p + i, j, :] for i in 1:K, j in 1:q]...)
Λ_chain = Chains(Λ_pos_sam);
#pΛ = plot(Λ_chain)
In [24]:
Λ
Out[24]:
5-element Array{Array{Float64,2},1}:
 [0.130044 -0.0289222 … 0.305525 0.264754; 0.0688177 -0.179345 … 0.127981 0.1679; … ; 0.0939935 0.234345 … 0.354274 0.384877; -0.171442 -0.232616 … -0.268721 0.369819]     
 [0.185065 -0.0558782 … 0.0867885 -0.55239; 0.775192 -0.126688 … 0.0377201 -0.494974; … ; -0.358786 -0.774431 … 0.658787 -0.76833; 0.762212 -0.771457 … -0.522754 -0.716634]
 [0.691138 0.022232 … -0.805842 -0.0223612; -0.25037 0.309294 … 0.311741 0.55864; … ; 0.227933 -0.337321 … 0.863367 -0.20838; 0.179381 0.239937 … 0.780241 -0.686083]       
 [0.551298 1.05886 … -0.386289 -1.04071; -0.403397 0.717535 … -1.08726 0.704522; … ; 0.395563 0.251704 … -0.603928 -0.66962; 0.385119 0.647842 … -1.1866 -1.46466]          
 [0.782915 0.473769 … 0.790897 -0.297221; -0.162902 1.29429 … 1.64004 -1.00375; … ; -1.31956 -0.271109 … 1.96226 1.76631; 1.15007 1.49836 … -1.51572 -1.7759]
In [25]:
ϕ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, K, 1);
ϕ_pos_sam[:, :, 1] = hcat([ϕ_sam[i, :] for i in 1:K]...);
ϕ_chain = Chains(ϕ_pos_sam);
pϕ = plot(ϕ_chain)
Out[25]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 4 6 8 Param1 Iteration Sample value 2 4 6 8 0.0 0.1 0.2 0.3 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Param2 Iteration Sample value 0 5 10 15 20 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 4 6 8 10 12 Param3 Iteration Sample value 4 6 8 10 12 0.00 0.05 0.10 0.15 0.20 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 5 10 15 20 25 Param4 Iteration Sample value 5 10 15 20 25 0.000 0.025 0.050 0.075 0.100 Param4 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 4 6 8 10 Param5 Iteration Sample value 2 4 6 8 10 0.0 0.1 0.2 0.3 Param5 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2.5 5.0 7.5 10.0 Param6 Iteration Sample value 2.5 5.0 7.5 10.0 12.5 0.00 0.05 0.10 0.15 Param6 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 4 6 8 10 12 14 Param7 Iteration Sample value 3 6 9 12 15 0.00 0.05 0.10 0.15 Param7 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 10 15 20 25 Param8 Iteration Sample value 5 10 15 20 25 30 0.00 0.02 0.04 0.06 0.08 Param8 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2.5 5.0 7.5 10.0 12.5 15.0 17.5 Param9 Iteration Sample value 0 5 10 15 20 0.00 0.05 0.10 0.15 0.20 Param9 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 5 10 15 20 25 Param10 Iteration Sample value 0 10 20 30 0.00 0.02 0.04 0.06 Param10 Sample value Density
In [26]:
Σ_pos_sam = Array{Float64, 3}(undef, N_sam + 1, q, 1);
Σ_pos_sam[:, :, 1] = hcat([Σ_sam[ind, :] for ind in 1:q]...);
Σ_chain = Chains(Σ_pos_sam);
pΣ = plot(Σ_chain)
Out[26]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 Param1 Iteration Sample value 0 5 10 15 20 0.0 0.5 1.0 1.5 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 Param2 Iteration Sample value 0 5 10 15 0.0 0.5 1.0 1.5 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 25 Param3 Iteration Sample value 0 5 10 15 20 25 0 1 2 3 4 5 Param3 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 10 20 30 Param4 Iteration Sample value 0 10 20 30 0.0 0.2 0.4 0.6 0.8 Param4 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 Param5 Iteration Sample value 0 5 10 15 20 0 1 2 3 4 5 Param5 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 Param6 Iteration Sample value 0 5 10 15 20 0.00 0.25 0.50 0.75 1.00 Param6 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 25 Param7 Iteration Sample value 0 10 20 30 0.0 0.2 0.4 0.6 Param7 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 20 Param8 Iteration Sample value 0 5 10 15 20 0 1 2 3 4 5 Param8 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 10 20 30 Param9 Iteration Sample value 0 10 20 30 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Param9 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 5 10 15 Param10 Iteration Sample value 0 5 10 15 0 1 2 3 Param10 Sample value Density
In [27]:
sqrt_Σ_diag.^2
Out[27]:
10-element Array{Float64,1}:
 0.5000000000000001 
 1.0                
 0.4                
 2.0000000000000004 
 0.29999999999999993
 2.5000000000000004 
 3.5                
 0.45               
 1.4999999999999998 
 0.5000000000000001
In [28]:
ω_incp_obs_pos_sam = Array{Float64, 3}(undef, n, q, N_sam);
lll = fill(1.0, (n, 1));
for i in 1:N_sam
    ω_incp_obs_pos_sam[:, :, i] = F_sam[:, :, i] * γ_sam[(p + 1):(p + K), :, i + 1] + 
                     lll * transpose(γ_sam[1, :, i + 1]);
end
truncindex = 1;#Integer(trunc(N_sam / 2));
ω_incp_pos_sam = Array{Float64, 3}(undef, N_sam - truncindex  + 1, 3, 1);
ω_incp_pos_sam[:, :, 1] = hcat(ω_incp_obs_pos_sam[1, 1, truncindex:N_sam], 
    ω_incp_obs_pos_sam[1, 2, truncindex:N_sam], ω_incp_obs_pos_sam[200, 1, truncindex:N_sam]);
ω_incp_chain = Chains(ω_incp_pos_sam);
pΣ = plot(ω_incp_chain)
Out[28]:
0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 2 4 6 8 Param1 Iteration Sample value 2 4 6 8 0.0 0.2 0.4 0.6 Param1 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 0 2 4 6 Param2 Iteration Sample value -2 0 2 4 6 0.0 0.1 0.2 0.3 0.4 Param2 Sample value Density 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 -2 -1 0 1 2 3 4 Param3 Iteration Sample value -2 0 2 4 0.0 0.2 0.4 0.6 Param3 Sample value Density
In [29]:
# check the variance covariance across latent process
cov_pos = Array{Float64, 3}(undef, q, q, N_sam);
for i in 1:N_sam
    cov_pos[:, :, i] = cov(F_sam[:, :, i] * γ_sam[(p + 1):(p + K), :, i + 1])
end   
# ω_incp_obs_pos_sam[:, :, i]
cov_pos_sam = Array{Float64, 3}(undef, N_sam, q * q, 1);
cov_pos_sam[:, :, 1] = hcat([cov_pos[i, j, :] for i in 1:q, j in 1:q]...);
cov_pos_chain = Chains(cov_pos_sam);
#pcov = plot(cov_pos_chain)
In [30]:
covω = cov(ω_ord[S, :])
Out[30]:
10×10 Array{Float64,2}:
  25.3808     7.60466     9.79245   …  -0.859226  -18.197     -11.1629  
   7.60466   18.2329      2.28601      -0.181906   -6.16284    -6.00531 
   9.79245    2.28601    29.5927        0.248116  -15.4448    -10.8923  
  -2.04999    4.72879    -6.63682      -4.7299      1.42522    -0.550013
  -1.25728    7.50492    -2.7964        2.56391    -1.73796    -0.341501
  12.1198     2.51267    11.7798    …   0.435206  -10.2852    -11.9809  
 -10.5902    -6.36495     5.03186      -0.673167   -0.655199   -0.213824
  -0.859226  -0.181906    0.248116     25.1696      9.2675      2.94299 
 -18.197     -6.16284   -15.4448        9.2675     34.803      12.8841  
 -11.1629    -6.00531   -10.8923        2.94299    12.8841     25.3644

Posterior Inference

In [31]:
N_Inf_burn = Integer(trunc(0.75 * N_sam));
ω_incp_obs_pos_qt = Array{Float64, 3}(undef, n, q, 3);
for j in 1:q
    for i in 1:n
        ω_incp_obs_pos_qt[i, j, :] = quantile(ω_incp_obs_pos_sam[i, j, N_Inf_burn:N_sam], [0.025, 0.5, 0.975])
    end
end
# count the covarage of 95% CI #
count_ω_incp = fill(0.0, q);
for j in 1:q
    for i in 1:n
        count_ω_incp[j] = count_ω_incp[j] + 
        ((ω_incp_obs_pos_qt[i, j, 1] < ω_incp_obs[S[i], j]) && 
            (ω_incp_obs_pos_qt[i, j, 3] > ω_incp_obs[S[i], j]))
    end
end
count_ω_incp
Out[31]:
10-element Array{Float64,1}:
 1073.0
 1163.0
  994.0
 1096.0
 1091.0
 1119.0
 1114.0
 1004.0
 1118.0
 1065.0
In [32]:
count_ω_incp ./ n
Out[32]:
10-element Array{Float64,1}:
 0.8941666666666667
 0.9691666666666666
 0.8283333333333334
 0.9133333333333333
 0.9091666666666667
 0.9325            
 0.9283333333333333
 0.8366666666666667
 0.9316666666666666
 0.8875
In [33]:
sum(count_ω_incp) / (q*n)
Out[33]:
0.9030833333333333
In [34]:
summary_table = Array{Float64, 2}(undef, 15, 5);
summary_table[1, :] = vcat(β[2, 1], mean(γ_sam[2, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[2, :] = vcat(β[3, 1], mean(γ_sam[3, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[3, :] = vcat(β[2, 2], mean(γ_sam[2, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[4, :] = vcat(β[3, 2], mean(γ_sam[3, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[5, :] = vcat(β[2, 3], mean(γ_sam[2, 3, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 3, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[6, :] = vcat(β[3, 3], mean(γ_sam[3, 3, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 3, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[7, :] = vcat(sqrt_Σ_diag[1]^2, mean(Σ_sam[1, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[8, :] = vcat(sqrt_Σ_diag[2]^2, mean(Σ_sam[2, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[9, :] = vcat(sqrt_Σ_diag[3]^2, mean(Σ_sam[3, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[3, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[10, :] = vcat(sqrt_Σ_diag[4]^2, mean(Σ_sam[4, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[4, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[11, :] = vcat(covω[1, 1], mean(cov_pos_sam[N_Inf_burn:N_sam, 1, 1]),
    quantile(cov_pos_sam[N_Inf_burn:N_sam, 1, 1], [0.5, 0.025, 0.975]));
summary_table[12, :] = vcat(covω[1, 2], mean(cov_pos_sam[N_Inf_burn:N_sam, 2, 1]),
    quantile(cov_pos_sam[N_Inf_burn:N_sam, 2, 1], [0.5, 0.025, 0.975]));
summary_table[13, :] = vcat(ϕ[1], mean(ϕ_sam[1, N_Inf_burn:N_sam]),
    quantile(ϕ_sam[1, N_Inf_burn:N_sam], [0.5, 0.025, 0.975]));
summary_table[14, :] = vcat(ϕ[2], mean(ϕ_sam[2, N_Inf_burn:N_sam]),
    quantile(ϕ_sam[2, N_Inf_burn:N_sam], [0.5, 0.025, 0.975]));
summary_table[15, :] = vcat(ϕ[3], mean(ϕ_sam[3, N_Inf_burn:N_sam]),
    quantile(ϕ_sam[3, N_Inf_burn:N_sam], [0.5, 0.025, 0.975]));
summary_table = round.(summary_table; digits = 3);
rnames = ["β[2, 1]", "β[2, 2]", "β[3, 1]", "β[3, 2]", "β[4, 1]", "β[4, 2]", "Σ[1, 1]", 
    "Σ[2, 2]", "Σ[3, 3]", "Σ[4, 4]", "cov(ω)[1, 1]", "cov(ω)[1, 2]", "ϕ1", "ϕ2", "ϕ3"];
summary_table = [rnames summary_table];
pretty_table(summary_table,  ["" "true" "mean" "median" "2.5%" "97.5%"], markdown)

|              |   true |    mean |  median |    2.5% |  97.5% |
|--------------|--------|---------|---------|---------|--------|
|      β[2, 1] |   -5.0 |   -5.14 |  -5.143 |  -5.563 | -4.714 |
|      β[2, 2] |    8.0 |   8.343 |   8.338 |   7.889 |  8.841 |
|      β[3, 1] |    2.0 |   2.166 |   2.164 |   1.696 |  2.637 |
|      β[3, 2] |    6.9 |   6.684 |   6.681 |   6.212 |  7.166 |
|      β[4, 1] |    3.0 |   2.856 |   2.855 |   2.464 |  3.239 |
|      β[4, 2] |  -12.0 | -12.153 | -12.152 | -12.571 | -11.72 |
|      Σ[1, 1] |    0.5 |   0.277 |   0.255 |   0.097 |  0.586 |
|      Σ[2, 2] |    1.0 |   1.243 |   1.225 |   0.812 |   1.78 |
|      Σ[3, 3] |    0.4 |   0.126 |   0.115 |   0.062 |  0.241 |
|      Σ[4, 4] |    2.0 |   1.566 |   1.523 |    0.81 |  2.458 |
| cov(ω)[1, 1] | 25.381 |  25.589 |  25.603 |  24.928 | 26.218 |
| cov(ω)[1, 2] |  7.605 |   7.538 |   7.537 |   7.022 |   8.06 |
|           ϕ1 | 16.792 |   4.956 |   4.921 |   3.222 |  6.418 |
|           ϕ2 | 20.078 |  13.949 |  14.216 |   8.636 | 17.334 |
|           ϕ3 | 14.988 |   8.036 |   8.316 |   4.613 | 10.538 |
In [35]:
summary_table = Array{Float64, 2}(undef, 20, 5);
summary_table[1, :] = vcat(β[2, 1], mean(γ_sam[2, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[2, :] = vcat(β[2, 2], mean(γ_sam[2, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[3, :] = vcat(β[2, 3], mean(γ_sam[2, 3, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 3, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[4, :] = vcat(β[2, 4], mean(γ_sam[2, 4, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 4, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[5, :] = vcat(β[2, 5], mean(γ_sam[2, 5, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 5, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[6, :] = vcat(β[2, 6], mean(γ_sam[2, 6, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 6, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[7, :] = vcat(β[2, 7], mean(γ_sam[2, 7, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 7, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[8, :] = vcat(β[2, 8], mean(γ_sam[2, 8, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 8, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[9, :] = vcat(β[2, 9], mean(γ_sam[2, 9, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 9, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[10, :] = vcat(β[2, 10], mean(γ_sam[2, 10, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[2, 10, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[11, :] = vcat(β[3, 1], mean(γ_sam[3, 1, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[12, :] = vcat(β[3, 2], mean(γ_sam[3, 2, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[13, :] = vcat(β[3, 3], mean(γ_sam[3, 3, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 3, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[14, :] = vcat(β[3, 4], mean(γ_sam[3, 4, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 4, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[15, :] = vcat(β[3, 5], mean(γ_sam[3, 5, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 5, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[16, :] = vcat(β[3, 6], mean(γ_sam[3, 6, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 6, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[17, :] = vcat(β[3, 7], mean(γ_sam[3, 7, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 7, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[18, :] = vcat(β[3, 8], mean(γ_sam[3, 8, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 8, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[19, :] = vcat(β[3, 9], mean(γ_sam[3, 9, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 9, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[20, :] = vcat(β[3, 10], mean(γ_sam[3, 10, N_Inf_burn:(N_sam + 1)]),
    quantile(γ_sam[3, 10, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table = round.(summary_table; digits = 3);
rnames = ["β[2, 1]", "β[2, 2]", "β[2, 3]", "β[2, 4]", "β[2, 5]", "β[2, 6]", "β[2, 7]", 
    "β[2, 8]", "β[2, 9]", "β[2, 10]", "β[3, 1]", "β[3, 2]", "β[3, 3]", "β[3, 4]", 
    "β[3, 5]", "β[3, 6]", "β[3, 7]", "β[3, 8]", "β[3, 9]", "β[3, 10]"];
summary_table = [rnames summary_table];
pretty_table(summary_table,  ["" "true" "mean" "median" "2.5%" "97.5%"], markdown)

|          |  true |    mean |  median |    2.5% |  97.5% |
|----------|-------|---------|---------|---------|--------|
|  β[2, 1] |  -5.0 |   -5.14 |  -5.143 |  -5.563 | -4.714 |
|  β[2, 2] |   2.0 |   2.166 |   2.164 |   1.696 |  2.637 |
|  β[2, 3] |   3.0 |   2.856 |   2.855 |   2.464 |  3.239 |
|  β[2, 4] |  -2.0 |  -1.347 |  -1.351 |  -1.934 | -0.751 |
|  β[2, 5] |  -6.0 |  -5.743 |  -5.749 |  -6.219 | -5.226 |
|  β[2, 6] |   4.0 |   4.156 |   4.152 |    3.61 |  4.711 |
|  β[2, 7] |   5.0 |   5.162 |   5.159 |   4.487 |  5.802 |
|  β[2, 8] |  -3.0 |  -2.924 |  -2.918 |  -3.361 | -2.512 |
|  β[2, 9] |   6.0 |   5.951 |    5.95 |   5.367 |  6.524 |
| β[2, 10] |  -4.0 |  -3.973 |   -3.97 |  -4.425 | -3.566 |
|  β[3, 1] |   8.0 |   8.343 |   8.338 |   7.889 |  8.841 |
|  β[3, 2] |   6.9 |   6.684 |   6.681 |   6.212 |  7.166 |
|  β[3, 3] | -12.0 | -12.153 | -12.152 | -12.571 | -11.72 |
|  β[3, 4] |   0.0 |   0.123 |   0.119 |   -0.46 |  0.714 |
|  β[3, 5] |  -4.0 |  -3.907 |  -3.917 |   -4.33 |  -3.42 |
|  β[3, 6] |   7.7 |   7.924 |   7.925 |   7.391 |  8.464 |
|  β[3, 7] |  -8.8 |  -8.791 |  -8.789 |  -9.458 | -8.114 |
|  β[3, 8] |   3.3 |    3.73 |   3.731 |   3.194 |  4.202 |
|  β[3, 9] |   6.6 |   6.651 |   6.651 |   6.128 |   7.18 |
| β[3, 10] |  -5.5 |  -5.301 |  -5.295 |  -5.766 |  -4.87 |
In [36]:
summary_table = Array{Float64, 2}(undef, 10, 5);
summary_table[1, :] = vcat(sqrt_Σ_diag[1]^2, mean(Σ_sam[1, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[1, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[2, :] = vcat(sqrt_Σ_diag[2]^2, mean(Σ_sam[2, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[2, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[3, :] = vcat(sqrt_Σ_diag[3]^2, mean(Σ_sam[3, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[3, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[4, :] = vcat(sqrt_Σ_diag[4]^2, mean(Σ_sam[4, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[4, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[5, :] = vcat(sqrt_Σ_diag[5]^2, mean(Σ_sam[5, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[5, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[6, :] = vcat(sqrt_Σ_diag[6]^2, mean(Σ_sam[6, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[6, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[7, :] = vcat(sqrt_Σ_diag[7]^2, mean(Σ_sam[7, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[7, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[8, :] = vcat(sqrt_Σ_diag[8]^2, mean(Σ_sam[8, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[8, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[9, :] = vcat(sqrt_Σ_diag[9]^2, mean(Σ_sam[9, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[9, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table[10, :] = vcat(sqrt_Σ_diag[10]^2, mean(Σ_sam[10, N_Inf_burn:(N_sam + 1)]),
    quantile(Σ_sam[10, N_Inf_burn:(N_sam + 1)], [0.5, 0.025, 0.975]));
summary_table = round.(summary_table; digits = 3);
rnames = ["Σ[1, 1]", "Σ[2, 2]", "Σ[3, 3]", "Σ[4, 4]", "Σ[5, 5]", "Σ[6, 6]", "Σ[7, 7]", 
    "Σ[8, 8]", "Σ[9, 9]", "Σ[10, 10]"];
summary_table = [rnames summary_table];
pretty_table(summary_table,  ["" "true" "mean" "median" "2.5%" "97.5%"], markdown)

|           | true |  mean | median |  2.5% | 97.5% |
|-----------|------|-------|--------|-------|-------|
|   Σ[1, 1] |  0.5 | 0.277 |  0.255 | 0.097 | 0.586 |
|   Σ[2, 2] |  1.0 | 1.243 |  1.225 | 0.812 |  1.78 |
|   Σ[3, 3] |  0.4 | 0.126 |  0.115 | 0.062 | 0.241 |
|   Σ[4, 4] |  2.0 | 1.566 |  1.523 |  0.81 | 2.458 |
|   Σ[5, 5] |  0.3 | 0.177 |  0.156 | 0.073 | 0.412 |
|   Σ[6, 6] |  2.5 | 1.869 |  1.839 |  1.27 | 2.676 |
|   Σ[7, 7] |  3.5 |  2.46 |  2.437 | 1.548 | 3.485 |
|   Σ[8, 8] | 0.45 | 0.183 |  0.173 | 0.084 | 0.344 |
|   Σ[9, 9] |  1.5 | 1.047 |  1.021 | 0.559 | 1.653 |
| Σ[10, 10] |  0.5 | 0.247 |  0.211 | 0.096 | 0.615 |
In [37]:
N_Inf_burn
Out[37]:
15000
In [38]:
N_Inf_burn = Integer(trunc(0.75 * N_sam));
count_Y_M = fill(0.0, q);
Y_m_pos_qt = [Array{Float64, 2}(undef, length(M_ind[ind]), 3) for ind in 1:q];
Y_m_pos_mean = [Array{Float64, 1}(undef, length(M_ind[ind])) for ind in 1:q]
for i in 1:q
    for j in 1:length(M_ind[i])
        Y_m_pos_qt[i][j, :] = quantile(Y_m_sam[i][j, N_Inf_burn:N_sam], [0.025, 0.5, 0.975]);
        Y_m_pos_mean[i][j] = mean(Y_m_sam[i][j, N_Inf_burn:N_sam])
    end
end
for i in 1:q
    for j in 1:length(M_ind[i])
        count_Y_M[i] = count_Y_M[i] + ((Y_m_pos_qt[i][j, 1] < Y_ord[M_ind[i][j], i]) && 
         (Y_m_pos_qt[i][j, 3] > Y_ord[M_ind[i][j], i]))
    end
end
count_Y_M
Out[38]:
10-element Array{Float64,1}:
 193.0
 192.0
 188.0
 193.0
 187.0
 183.0
 192.0
 193.0
 191.0
 185.0
In [39]:
count_Y_M ./ 200
Out[39]:
10-element Array{Float64,1}:
 0.965
 0.96 
 0.94 
 0.965
 0.935
 0.915
 0.96 
 0.965
 0.955
 0.925
In [40]:
sum(count_Y_M)/ (q*200)
Out[40]:
0.9485
In [41]:
SPE = fill(0.0, q)
for i in 1:q
    for j in 1:length(M_ind[i])
        SPE[i] = SPE[i] + (Y_m_pos_mean[i][j] -  Y_ord[M_ind[i][j], i])^2
    end
end
sqrt.(SPE ./ 200 )
Out[41]:
10-element Array{Float64,1}:
 1.9016701308984667
 2.039719798001778 
 2.06981291517544  
 2.911048771709707 
 2.26006831205493  
 2.6935341074755788
 3.0058450676574457
 2.268021986574122 
 2.416107961362959 
 2.4909946662305975
In [42]:
sqrt(sum(SPE) / (200*q))
Out[42]:
2.4314570070226518
In [43]:
N_Inf_burn = Integer(trunc(0.75 * N_sam));
count_cov = 0.0
for i in 1:q
    for j in i:q
        count_cov = count_cov + ((quantile(cov_pos[i, j, N_Inf_burn:N_sam], [0.025])[1] < (covω[i, j])) 
         && (quantile(cov_pos[i, j, N_Inf_burn:N_sam], [0.975])[1] > (covω[i, j])))
    end
end
count_cov
Out[43]:
47.0
In [44]:
#@save "../results/ω_incp_obs_pos_qt_BSLMC.jld"  ω_incp_obs_pos_qt